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We investigate to what extend the replica-exchange Monte Carlo method is able to equilibrate 
a simple liquid in its supercooled state. We find that this method does indeed allow to generate 
accurately the canonical distribution function even at low temperatures and that its efficiency is 
about 10-100 times higher than the usual canonical molecular dynamics simulation. 
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If a liquid is cooled to a temperature close to its glass 
transition temperature, its dynamical properties show a 
drastic slowing-down. At the same time, a crossover from 
highly unharmonic liquid- like behavior to harmonic solid- 
like behavior is expected in its static (thermodynamic) 
properties at a certain temperature Tk, the Kauzmann 
temperature Very recently the value of Tk of sim- 
ple model liquids have been determined analytically |^] 
and numerically |^,^ and some possibilities of a ther- 
modynamic glass transition at Tk have been discussed. 
Although the values of Tk obtained with the different 
methods are consistent with each other, it was necessary 
for the numerical calculations of Tk to extrapolate high 
temperature data (T > 0.45) of the liquid and disordered 
solid branches of the configurational entropy S{T) down 
to significantly lower temperatures {Tk — 0.3). With a 
guide of an analytic prediction for liquids, S{T) ~ 
0, and for harmonic solids, S{T) ~ logT, a crossing of 
the two branches has been found and used to calculate 
Tk- However, to make those observations more reliable, 
very accurate calculations of thermodynamics properties 
are necessary in the deeply supercooled regime, which is 
difRculat since the typical relaxation times of the system 
are large. 

In recent years, several efficient simulation algorithms 
have been developed to generate canonical distributions 
also for complex systems. Examples are the multi- 
canonical 1^,0, the simulated tempering and the 
replica-exchange (RX) methods. Although these 
methods were originally developed for Ising-type spin sys- 
tems, their applications to any off-lattice model by use 
of Monte Carlo or molecular dynamics simulations are 
rather straightforward [p^-p^. However, it has been 
found that the application of some of these algorithms 
to supercooled liquids or structural glasses is of only lim- 
ited use p^. The main motivation of the present paper 
is to test the efficiency of the RX method, which seems 
to be in many cases the most efficient algorithm, to the 
case of highly supercooled liquids . 

The system we study is a two-component {AB) 
Lennard-Jones mixture, which is a well characerized 



model system for supercooled simple liquids. The to- 
tal number of particles is = 1000, and they inter- 
act via the (truncated and shifted) potential 4'ap{''^ij) = 
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where r,,- is the distance 
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between particles i and j, and the interaction parame- 
ters are a,/3 S A,B, caa = 1, ^ab = 1-5, tsB = 0.5, 
gaa = 1, <^AB = 0.8, and a be — 0.88. Other simulation 
parameters and units are identical as in | [T^ . The time 
step At for numerical integration is 0.018. 

The algorithm of our replica-exchange molecular dy- 
namics (RXMD) simulation is essentially equivalent to 
that of Ref. ]l5[|, and therefore we summarize our sim- 
ulation procedure only briefly, (i) We construct a sys- 
tem consisting of AI noninteracting subsystems (repli- 
cas), each composed of N particles, with a set of arbi- 
trary particle configurations {q^, • • ■ , ^^"^ momenta 
{Pi, ■ ■ ■ ,Pm}- The Hamiltonian of the m-th subsystem 
is given by 



) K{p„J + A„,E{q 
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where K is the kinetic energy, E is the potential energy, 
and Am e {Ai, • • • , Am} is a parameter to scale the po- 
tential, (ii) A MD simulation is done for the total sys- 
tem, whose Hamiltonian is given by 7i = J2m=i ^"^^ 
at a constant temperature T = using the con- 

straint method ||l^. Step (ii) generates a canonical 
distribution P(qi,---,qj,j;/3o) = H^„LiP(q,„; A^/Jq) oc 
exp[-l3Qj2m=i^'mE{q^)] in configuration space [|o). 
(iii) At each time interval AtB.x, the exchange of the po- 
tential scaling parameter of the m-th and 7i-th subsystem 
are considered, while {q^, • • • , q^} ^'^d {pi, ■ ■ ■ ,Pm} are 
unchanged. The acceptance of the exchange is decided in 
such a way that it takes care of the condition of detailed 
balance. Here we use the Metropolis scheme, and thus 
the acceptance ratio is given by 



1, 

exp(-A„i,„), 



A™,„ < 
A™,„ > 0, 



(2) 



where A,„,„ = /9o(A„ - A„)(£:(q„J - E{qJ). (iv) Re- 
peat steps (ii) and (iii) for a sufficient long time. This 
scheme leads to canonical distribution functions P{E; (ii) 
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at a set of inverse temperatures f3i = Ai/3o- To make a 
measurement at an inverse temperature f3i one has to av- 
erage over all those subsystems {i £ 1, • • ■ , M) for which 
we have (temporarily) pi = Ai/3o. Usual canonical molec- 
ular dynamics (CMD) simulations are realized if we skip 
step (iii). 

In the present simulation, we take M = 16, Po = 
0.45^^, Xi — 1 — 0.0367(i — 1) and thus cover a tempera- 
ture range 0.45 <T<1. Exchange events are examined 
only between subsystems that have scaling parameters 
Xi and Xij^i that are nearest neighbors; the events with 
i = 1, 3, 5, • • • or i = 2, 4, 6, • • • are repeated alternatively 
every Atjix intervals. We find that the highest average 
acceptance ratio for this type of move is 0.186 for the ex- 
change of Ai and A2, and the lowest is 0.027 for A15 and 
A16. Although these values can be made more similar by 
optimizing the different gaps between Ai and A^+i for a 
fixed choice of Ai and Xm, only small improvements were 
obtained by such a simple optimization in our case. We 
also note that the choice of Atax strongly affects the ef- 
ficiency of the RX method; At fix should be neither too 
smah or too large We used AtRx = 10^ At, a time 
which is a bit larger than the one needed for a particle to 
do one oscillation in its cage, and data are accumulated 
for 0<t<5xlO^At after having equilibrated the sys- 
tem for the same amount of time. At the beginning of 
the production run, the subsystems were renumbered so 
that at t = we had for all m A^ = A™. 

In Fig. 1(a), we show the time evolution of the subsys- 
tems in temperature space. One can see that the subsys- 
tems starting from the lowest (m = 1) and the highest 
(m = 16) temperature explore both the whole tempera- 
ture space from i ~ 1 to 16. Fig. 1(b) presents the mean 
squared displacements (MSD) 



Ai?'(t)-|9„(t)-g„(0)|VA^ 
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for the RXMD (with m = 1) and for the CMD performed 
at T = (Ai/3o)^^ ~ 0.45. From this figure we recognize 
that, due to the temperature variation in the RXMD 
method, the system moves very efficiently in configura- 
tion space, while in the CMD the system is trapped in a 
single metastable configuration for a very long time. If 
one uses the MSD to calculate an effective diffusion con- 
stant, one finds that this quantity is around 100 times 
larger in the case of the RXMD than in the CMD case, 
thus demonstrating the efficiency of the former method. 

Fig. 2 shows the canonical distribution function of the 
total potential energy at the different temperatures. 



P,{E)=P{E-XM, 
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obtained by a single RXMD simulation. For adjacent 
temperatures the corresponding distribution functions 
should have enough overlap to obtain a reasonable ex- 
change probabilities and hence can be used to optimize 



the efficiency of the algorithm. Further use of these dis- 
tribution functions can be made by using them to check 
whether or not one has indeed equilibrated the system. 
Using the reweighting procedure ]2^ , it is in principle 
possible to calculate the canonical distribution functions 



P^{E■,X,Po) 



P,{E) exp[{X, ~Xj)PoE] 
J dE' P,{E')exp[{X,-X,)poE'] 
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at a new temperature Tj = (XjPo)^^ from any Pi(E). 
Note that in equilibrium the left hand side should be 
independent of i to within the accuracy of the data. 

In Fig. 3 we plot different Pi(E; Xi^Po), using as input 
the distributions P{E; XiPo) for 1 < i < 8, obtained from 
RXMD (a) and CMD (b) simulations. (Both simulations 
extended over 8.7 x 10* time units.) Wc see that in the 
case of the RXMD the different distributions Pi fall nicely 
on top of each other in the whole energy range, thus giv- 
ing evidence that the system is indeed in equilibrium. In 
contrast to this, the different distributions of the CMD, 
Fig 3b, do not superimpose at low energies (=low tem- 
peratures), thus demonstrating the lack of equilibration. 
This can be seen more clealy by comparing Fig. 3(c) and 
(d), where Pi{E; XiPo) is plotted. 

Fig. 4 shows the temperature dependence of the po- 
tential energy E{T) obtained from RXMD simulations 



E{T,) = j dE' P{E';X,po)E' 
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For the sake of comparison we have also included in 
this plot data from CMD with the same length of the 
production run as well as data from CMD simulations 
which were significantly longer (about one order of mag- 
nitude) ^. The solid line is a fit to the RXMD results 
with the function E{T) ^ Eq + AT° '', a functional form 
suggested by analytical calculations Q . One can see that 
RXMD and CMD results coincide at higher tempera- 
tures, but deviations become significant at low tempera- 
tures (see Inset). Furthermore, we see that the present 
RXMD results agree well with CMD data of the longer 
simulations. 

As a final check to see whether the RXMD is indeed 
able to equilibrate the system also at low temperatures, 
we have calculated the temperature dependence of the 
(constant volume) heat capacity (T) via the two routes 



Cy{T) = dE{T)/dT 

= {{E') - {E^))IT\ 



(7) 
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and plot the results in Fig. 5. Again we see that within 
the accuracy of our data the two expressions give the 
same answer, thus giving evidence that the system is in- 
deed in equilibrium. 

Summary: We have done replica-exchange molecular 
dynamics and canonical molecular dynamics simulations 
for a binary Lennard-Jones mixture in order to check 
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the efficiency of the rephca-exchange method for a struc- 
tural glass former in the strongly supercooled regime. We 
find that at low temperatures the RXMD is indeed sig- 
nificantly more efficient than the CMD, in that the ef- 
fective diffusion constant of the particles is around 100 
times larger in the RXMD. However, accurate simula- 
tions are still difficult for T < 0.45 even with RXMD. 
Finding an optimal choice of M, {Ai, • • • , Am}, and AIrx 
may be important in order to allow simulations also for 
T < 0.45 within reasonable computation times. Further- 
more it might be that the efficiency of RXMD improves 
even more if one uses it below the critical temperature 
of mode-coupling theory , since there is evidence that 
below this temperature the nature of the energy land- 
scape is not changing anymore . 
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FIG. 1. (a) Typical walks of the subsystems in temperature space, (b) Time dependence of the mean squared displacement. 
The solid line shows Ai?^(t) from RXMD for a subsystem which at f = was at T = 0.45 {i — 1), and the dashed line is 
Aii^(t) from CMD at T = 0.45. The two curves have been calculated by starting from the same initial configuration. 




FIG. 2. The canonical distribution function Pi{E) at various temperatures Ti {1 < i < 16) obtained by a single RXMD 
simulation. Here, Ti = 0.45 and Tie = 1.0. 
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FIG. 3. The canonical distribution function at T = (A4/3o)~^ = 0.506 by reweighting Pr{E) for 1 < i < 8 obtained by RXMD 
(a) and standard CMD (b) simulations. Numbers in the parentheses present temperatures at which simulations were done. 
The same function at T = (Ai/3o)~^ = 0.45 obtained by RXMD (c) and CMD (d). Note that in both simulations the length of 
the runs is the same (8.7 x 10* time units). 
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FIG. 4. Temperature dependence of the potential energy E{T) obtained via RXMD(o) and CMD(+) of runs with the same 
length. * presents values from much longer CMD runs. The solid line is the best fit to the RXMD data with a fit function 
E = Eo + AT^ '^, where Eo = -8.656 and A = 2.639 are fit parameters. 
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FIG. 5. Temperature dependence of the heat capacity C„(r) obtained via RXMD. □ presents data from C„ = dE{T)/dT, 
and ■ presents data from Cv = {{E^) — {S)^)/r^. The solid line is the result of a fit C„ = 0.6AT~°'*, with the same value of 
A as in Fig. 4. 



